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(54) Method for determining server staffing 

(57) A method determines the number of servers as 
a function of time required for a finite server queueing 
system based on a projected load. The number of serv- 
ers is chosen subject to the constraint that the probabil- 
ity of delay before beginning service does not exceed a 
target probability at all times. The finite server queueing 
system is first advantageously modeled as an infinite- 
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server system and the mean and variance of the number 
of busy servers at time / is determined and a distribution 
is approximated. For any time, the number of servers is 
chosen to be the least number so that the probability 
that these are all busy is less than the target probability. 
As an optional refinement, the infinite server target tail 
probability is then related to the actual delay probability 
in the finite server system. 
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Description 
Technical Field 

5 The invention relates to the field of queueing or service systems. 

Background 

Queueing systems are systems in which customers arrive, wait their turn to be served, are served by a server, 
w and then depart. Supermarket checkout lines, telephone operator service centers, bank teller windows and doctor's 
waiting rooms are examples of queueing systems. Queueing systems are advantageously characterized by five com- 
ponents: 

1. the customer arrival process (formally described as a stochastic process, e.g. a random counting process that 
15 describes the probability that the number of customers arriving in any time interval is a value k), 

2. the service time process (more formally described by a probability cumulative distribution function for each time 
t which indicates the probability that the time required to serve a customer arriving at time t is less than or equal 
to a value y), 

3. the number of servers in the queueing system, 

20 4. the queueing discipline, which describes the order in which customers are taken from the queue (e.g. first come, 

first served or sickest patient first), and 

5. the amount of buffer space in the queue, which describes how long the queue may be before customers are 
"lost" or turned away from the system. 

25 See, Andrew S. Tancnbaum, Computer Networks, 2 cd., Prentice Hall, Englowood Cliffs, New Jersey. 1989; R.W. Hall, 
Queueing Methods for Services and Manufacturing, Prentice Hall, Englewood Cliffs, New Jersey. 1991. 

One important consideration in queueing systems is the determination of how many servers are needed as a 
function of time. In particular, it is often the case that servers are assigned as function of time in response to anticipated 
changing loads on the system (i.e., to projected changes in the customer arrival process or to projected changes in 

30 the service time required by customers). For example, the problem may be toassign a time-varying number of operators 
at a telephone office at the beginning of each day where the number of servers is allowed to change every hour or 
every quarter hour. This general problem is often referred to as the server staffing problem. 

It is advantageous to analyze the server staffing problem in queueing systems by partially characterizing compo- 
nents of a queueing system by one or more functions. In particular, the arrival process advantageously is at least 

35 partially characterized by an arrival rate, X(t) (i.e. , \(t)&t describes the expected number of arrivals in the time interval 
(t, t+At) for small Ay. If the arrival process is a Poisson process (useful in describing many chance phenomena), the 
arrival process is completely characterized by k(t). See, e.g., A. Papoulis, Probability, Random Variables, and Sto- 
chastic Processes, 2 ed., McGraw-Hill, New York, 1984. Otherwise, X(t) only partially characterizes the arrival process, 
and other functions, such as those describing variability about the rate X(t) need to be specified. 

40 The problem of server staffing based on projected loads is a difficult problem if the arrival rates and/or service time 

processes of the queueing system are not constant (e.g. when X(t) is not a constant for all t). If the arrival rates and/ 
or service time processes of the queueing system vary with time, then the queueing system is said to be nonstationary. 
Although the server staffing problem with a stationary arrival rate and stationary service time process has been analyzed 
extensively, if the queueing system is nonstationary, then it is difficult to determine the number of servers as a function 

45 of time, s(t), and approximations are typically used See, e.g., Hall, supra. 

One way to address the nonstationary server staffing problem is to use a steady state analysis. In such an analysis 
either 1 ) the arrival rate is assumed to change sufficiently slowly relative to the service times so that steady state 
formulas can be applied at specific points in time using the arrival rates prevailing at those times or 2) the arrival rate 
is assumed to be some long term average arrival rate and steady-state formulas are applied using this long term 

so average arrival rate. See, e.g., Hall, supra. 

Another way to address the nonstationary server staffing problem is to recognize that the probability distribution 
of the number of customers present at any given time is advantageously approximated as a normal distribution, which 
normal distribution is specified by mean and variance parameters. The mean and variance of the approximated normal 
distribution have been determined by various ad hoc methods, e.g., by refining to the steady state analysis. See, e.g., 

55 G.M. Thompson, "Accounting for the Multi- Period Impact of Service When Determining Employee Requirements for 
Labor Scheduling," J Opns. Mgmt, vol. 11, pp. 269-288 (1993). Given the probability distribution of the number of 
customers in the queueing system at time t it is possible to select the smallest possible value j for the number of 
servers such that the probability that all j servers are busy at time t is less than or equal to some target value y The 
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value y determines the probability that a customer arriving at time t will encounter a delay before beginning to receive 
service If p n is the probability that the number of customers present is a value n at time /. then / is selected such that 



However, the techniques described above for allocating servers in nonstationary queueing systems typically are 
w not sufficiently accurate to determine server staffing as a function of time. In particular the techniques typically over- 
allocate servers at some times (thereby increasing system costs since not all servers will be fully utilized) and/or under- 
allocate servers at other times (thereby causing congestion, i.e. long waiting times for customers in the queue). Hence 
there is a need for better techniques for allocating servers in nonstationary queueing systems. 

is Summary 

In accordance with the present invention it is recognized that the number of servers as a function of time in a 
nonstationary finite queueing system may be allocated based on a probability distribution of the number of busy servers 
in a nonstationary queueing system with an infinite number of servers. In particular, the probability distribution is char- 
20 acterized by a mean function and a variance function, wherein the mean and variance functions are expressed in terms 
of input parameters that characterize the finite server queueing system. 

Brief Description of the Drawings 

25 FIG. 1 illustrates a system in which the inventive method may be used. 

FIG. 2 illustrates determination of server levels based on a probability distribution. 
FIG. 3 illustrates steps in the inventive method. 

Detailed Description 



The detailed description is divided into three parts. In Part I an overview of the inventive method is presented. In 
particular, a system in which the inventive method may be used is described and the steps in the inventive method are 
explained. Part I also defines terms useful in characterizing the operation of queueing systems. Part II, describes 
techniques for implementing the inventive method and discusses the inventive method for queueing systems having 
35 specific characteristics. Part it! discusses other applications of the inventive method. 



FIG. 1 illustrates a queueing system in which the inventive method for allocating servers may be used. Queueing 



40 system 100 comprises queue 105 and service portion 120. Customers 101 -/'arrive and take a place in queue 105. J n 
this exemplary system queue 105 is advantageously assumed to be infinitely long (so customers are never turned 
away or lost), and the service discipline is advantageously assumed to be first come-first served (FCFS). The arrival 
process of customers to the queue is advantageously characterized by one or more functions. For example, if the 
arrival process is Poisson, a single function, the arrival rate function X(t) } can be used to characterize the arrival process 

45 exactly (then X(t)At describes the expected number of arrivals in the time interval (t, t+At) for small At). In other cases, 
the arrival process may be partially described by a time-dependent arrival rate function, X(t), and optionally by a time 
dependent variability function, c a 2 (t),which describes the variability in the (random) arrival process about its determin- 
istic rate \(t). With the functions and c g 2 (t), two kinds of variability are captured. One kind of variability is represented 
by changes in the rate k(t) over time f. A second kind of variability (uncertainty or randomness) is captured by the 

so function c a 2 (t). 

Referring again to FIG. 1, a customer exits the queue and enters service area 112-k in service portion 120 when 
server 115-fr is available. For example in FIG. 1, customer 101-1 will exit queue 105 and enter service area 112-3 
because server 115-3 is available after having finished providing service to customer 101-6 who has exited queue 
system 100. 

55 The service times for each customer are advantageously treated as mutually independent and as independent of 

the arrival process. A cumulative distribution function, G ( (x), is used to described the probability that the service time 
will be less than or equal to x given that a customer arrives at time t. As shown below, G t (x), if not known exactly, may 
advantageously be approximated based on its mean value. 
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The inventive method determines the required number of servers as a function of time, s(t), required to meet some 
criterion based on parameters characterizing the arrival and service processes, eg \(t), c a 2 (t)nx\6 G t (x) In particular, 
s(t) is advantageously determined so that the probability of delay (i.e., the probability that a customer arriving at time 
f would have to wait before beginning to receive service) is approximately equal to a target probability at all times. The 
s inventive method advantageously assumes that an infinite server queue is available, and the mean and variance of a 
random variable, QJt), describing the number of busy, servers at time / in that infinite server queue is determined. (A 
random variable is a quantity that takes on a value in a range of values with a certain probability, i.e., that has a 
probability cumulative distribution function). In particular, the inventive method advantageously assumes initially that 
an infinite number of servers are available in the system which is essentially an averaging procedure that produces an 
10 effective arrival rate at all times (when scaled by the mean service time). A probability cumulative distribution function 
for QJt) is then used to determine the number of servers as a function of time. In particular, the number of servers is 
advantageously chosen to be the smallest possible value that the probability that all these chosen servers would be 
busy at time t is less than or equal to a target value. 

FIG. 2 illustrates this process for an illustrative probability cumulative distribution function at a specific time t r 
15 Suppose the target probability is 0.1 , i.e., there is a 10% chance that a customer will have to wait or alternatively, a 
90% chance that customer will not have to wait. The idea to select a number of servers at time t r s(t r ), so that the value 
of the probability cumulative density for that selected number of servers is just greater than 0.9, i.e. the probability of 
not having to wait. Thus, referring to FIG. 2, the number of servers at time t r should be selected to be eleven. If the 
number of servers is jaken to be eleven, then the likelihood that more than eleven servers will be busy, i.e. that a 
20 customer will have to wail because the finite server system has only eleven servers, is slightly less than 0.1 . 

FIG. 3 illustrates steps in the inventive method for allocating servers in a queueing system, and each step is 
described in detail below. In step 310, parameters characterizing the arrival process and service time process are used 
to determine the mean of the number of busy servers at time f assuming an infinite server queue. For the infinite server 
queue, assumed to start empty at some time prior to /, the mean number of busy servers as a function of time is (exactly) 



where GP(t)= 1 - GJt)s P(S(u) > t) with S(u) being the service time of an arrival at time u and X(u) is the arrival rate 
at time u. Given only a mean y(t) for G/x| G/xj can be approximated by an exponential distribution with mean y(t); 
then G/x) is treated as (1~exp('X/y(t)). 

Equation (1) is known to be valid for the case of nonstationary Poisson arrivals, but equation (1) remains valid 

3S without the Poisson assumption, provided that the arrival-rate function X(t) is well defined. See, S.G. Eick, W. A. Massey 
and W. Whitt, "The Physics of the M/G/oo Queue, - Open Res., vol. 41 , pp. 731 -742 (1 993); W. A. Massey and W. Whitt, 
"Networks of Infinite Server Queues with Nonstationary Poisson I nput," Queueing Systems, vol. 13, pp. 183-250(1993). 
Thus, m(t) may be determined using equation (1 ), either analytically or numerically, depending on the nature of the 
functions X and G t Simplifying assumptions and/or approximations that may be used to facilitate the determination of 

40 m(t) are discussed in Part II below. 

In step 330 a determination is made of the variance of QJt). Note that if the arrival process is assumed to be 
Poisson, the number of busy servers at time t in the infinite server model has a Poisson distribution with the mean m 
(t) as in equation (1 ), and computations for step 330 are not necessary. Hence, as is the case with Poisson distributions, 
the variance v(t) equals the mean m(t). More generally, m(t) is regarded as the default value for v(t). If the arrival 

45 process is not Poisson, the variance vft) of QJt) is advantageously determined based on W. Whitt, "Understanding 
the Efficiency of Multi-Server Service Systems," Management Sci., vol. 38, pp. 708-723 (1992). In particular, the fol- 
lowing approximation is used: 



25 
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v(t) = z(t)m(t), 
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for some suitable r\ > 0 and where the expected value operator; in (4) Var is the variance, fl.e. VarX=E[X 2 l -(E 
[XJ) 2 ) and 4(7; counts the number of arrivals before time f. The function z(t) in (3) is a time-dependent generalization 
of the heavy-traffic peakedness. The peakedness is the ratio of the variance to the mean of the number of busy servers 

w in a stationary infinite server queueing system. The heavy-traffic approximation is the limiting value as the arrival rate 
is made to increase to infinity See, W. Whitt, "Heavy -Traffic Approximations for Service Systems with Blocking," AT& T 
Bell Lab. Tech. J., vol. 63, pp. 689-708 (1 984). The function c\(X) in (4) is a time-dependent generalization of the.index 
of dispersion for counts (the variance divided by the mean of the number of arrivals in an interval). See, K.W. Fendick 
and W. Whitt, "Measurementx and Approximations to Describe the Offered Traffic and Predict the Average Workload 

is in a Single-Server Queue,' Proc. IEEE, vol. 77, pp. 171-194 (1989). 

In step 340 once the mean and variance are found, an approximation for the distribution of the number of busy 
servers at time t, O^ft), is determined. Given the approximate distribution, the number of servers as a function of time, 
s(t), is selected in step 350 so that 



P(QJty> s(t)) <a (5a) 



2S 



and 

P(OJt) > (s(t)-1))>a (5b) 



for some prescribed target probability a. Equation (5a) represents a determination of a value of s(l) such that the 
probability of encountering a busy server at time f is less than a while equation (5b) represents a determination that 

30 one less server than the value of s(t) would render the probability of encountering a busy server at time f greater than a. 

Tho invontivo method rocognizes that a normal approximation may advantageously be used with the rcsuft that 
computations for (5a) and 5(b) are simplified. Let N(m, a 2 ) denote a normal random variable with mean mand variance 
o 2 Then the distribution of QJt) may be approximated by the distribution of N (m(t) f v(t)) for a suitable mean function 
m(t) and variance function v(t) as described above. Hence, [Q*>(t) - m(t)]t Jv(t) is distributed approximately as N(0, 1), 

35 so that the number of servers is taken to be 

s(r; = r m(t)+z a Jv(T) (6) 

40 where r xi is the least integer greater than xand z a satisfies 

P(N(0,1>zJ=a (7) 

45 Equations (6) and (7) simply define the number of servers necessary so that the area under the normal distribution for 
the number of servers is greater than the selected number at time t is less than the target probability value. It is also 
advantageous to incorporate a refinement to account for the discreteness {s(t) should be an integer) and thus 

so s(t) = r m(t) + 0.5 + z a JvU) (8) 

is advantageously used. 

Finally, in step 360 refinements to the function s(t) may optionally be made. The first refinement is based on known 
differences between finite-server and infinite server systems and the second refinement is based on a modified offered 
55 load approximation as described in detail below. 

Thus, the inventive method recognizes that the number of servers as a function of time may be allocated, subject 
to a criterion, based on a probability distribution of the number of busy servers, wherein the probability distribution is 
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characterized by a time-dependent mean function and a time dependent variance function, wherein the mean and 
variance functions are derived from a nonsfationary infinite server queueing system.. 

II Implementation of the Inventive Method 

5 

This section is divided into three parts The first part describes simplifying approximations and/or assumptions that 
may be used in step 310 of FIG. 3 to determine the mean, m(t), of OJt), the number of busy servers at time t in the 
infinite server system. The second part discusses further considerations in the determination of the variance, v(t), of 
O.Jt). The final part describes in detail refinements that may be applied in determining s(t) 

w 

A. Considerations in determining m(t) 

The inventive method (as illustrated in step 310 of FIG. 3) requires a determination of the mean, m(t), of QJt), the 
number of busy servers given an infinite number of servers at time t. Advantageously, equation (1 ) may be used, but 
'5 simplifying assumptions and/or. approximations can facilitate the determination. For example, if the service-time dis- 
tribution G,does not depend on t, then equation (1 ) can be rewritten as 



m(/) = E[ \X(u)du] = E{X(t - S t )]£[S], (9) 



where S is a generic service time and S e is a generic service-time stationary-excess random variable, i.e., 

25 

P(S< Z [P(S > u)du (10) 



30 

see Eick et al., supra. Note from formulas (9) and (10) that the infinite server system mean m(t) can be regarded as a 
smoothing of X(t) multiplied by the mean service time. Thus, the inventive method advantageously assumes initially 
that an infinite number of servers are available in the system which is essentially an averaging procedure that produces 
an effective arrival rate at all times (when scaled by the mean service time). 

3S in the case where service times are time dependent, an approximation for time-dependent service times can be 

developed, assuming that G/x) changes relatively slowly compared to X(t). There are many common situations in 
which this is the case. For example, measurements indicate that the holding-time distribution of telephone calls changes 
throughout the day, having a bigger mean in the evening than in the day, but this rate of change is slow compared to 
the rate of change of the arrival rate X(t). Assuming that the service-time distribution does indeed change relatively 

40 slowly, the following approximation to equation (9) based on equation (10) is introduced, namely, 



m(t)*ElMt-S 0 (t))]E[S(t)l (11) 
45 where S(t) has cumulative distribution function G t and S e (t) has cumulative distribution function 

I 



so 



p ( s AO^x) = 1 ^[P(S(t)>u)du, x>0. (12) 
E[S{t)] * 



Equation (1) above assumes that if the service-time distribution is associated with customers and is determined 
upon arrival. If instead each server is thought of as completing service at rate \i(t) at time t, then S(t) can be approximated 
by an exponential random variable with mean 1/ \x(t). If, however, the service rate changes relatively quickly, then the 
55 exact (time-dependent exponential) cumulative distribution function could be used, i.e., G t (x)=1-& t (x), where 
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G:(x) = e '^, *>0. (13) 

s To apply the inventive method staffing procedure described in the overview above, the mean m(t) as a function oft 

for all times t of interest needs to be computed. The formulas (6)-( 1 1 ) can be used directly, but they are not especially 
convenient for this purpose. However, if further assumptions are made, then very convenient procedures for computing 
m(t) can be obtained. In particular, efficient procedures can be obtained if simplifying assumptions about the service- 
time distribution are made. 

w First, if the service-time distribution is exponential with constant rate \i, then m(t) can be obtained simply by solving 

the ordinary differential equation (ODE) 

m7y=X(-y-mWn, (14) 

15 

where m(0) is given, e.g., m(0) = 0 it the queue is empty at time 0; see Eick et al., supra. Note that equation (14) 
extends to time-dependent service rates. More generally, if there is a fixed phase-type service-time distribution with k 
phases, then the mean m(t) can be obtained by solving a system of /cODEs. See, J.L. Davis, W.A. Massey and W. 
Whitt, "Sensitivity to the Service-Time Distribution in the Nonstationary Erlang Loss Model," Management Sci., vol. 42 
20 pp. 1107-1116, (1995); W.A. Massey and W. Whilt, "An Analysis of the Modified Offered Load Approximation for the 
Nonstationary Erlang Loss Model," Ana Appi Pcob., vol. 4, pp. 1145-1160 (1994). For background on phase-type 
distributions, see M.F. Neuts, Matrix-Geometric Solutions in Stochastic Models, The Johns Hopkins University Press, 
Baltimore, 1981. 

Second, instead of making special assumptions about the service-time distribution, it may be advantageous to 
2S assume a special structure for thc*arrival rate function X(t). Explicit formulas for m(t) when X(t) is sinusoidal arc given 
in S.G. Eick, W.A. Massey and W. Whitt, "M/G/<~ Queues with Sinusoidal Arrival Rates," Management Sci., vol. 39, 
pp. 241-252 (1993). Explicit formulas for m(t) when X(t) is a polynomial or a step function are given in S.G. Eick, W.A. 
Massey and W. Whitt, "The Physics of the M/G/~ Queue," Oper. Res., vol. 41 . pp. 731 -742 (1993). 

However, in practice X(t) typically will not have a well-defined mathematical structure. A convenient technique for 
30 determining m(t) when X(t) is a general function estimated from data is as follows. General arrival rate functions arising 
in practice are treated by making polynomial (e.g., linear or quadratic) approximations to the arrival-rate function X(t) 
over subintervals. For example, mftlcan be computed over the interval [(/c-1 )c. kc) by using a polynomial approximation 
to X(t) over the interval [(k-j)c, kc), where c is an appropriate constant such as 5 mean service times and j is perhaps 
about 8. The polynomial approximation for X(t) can be fit to data using least squares or iterated weighted least squares. 
35 (Iterated weighted least squares yields the maximum likelihood estimator, but ordinary least squares usually performs 
almost as well.) 

The ODEs provide a simple effective way to calculate m(t) for general X(t). Nevertheless, in some situations it may 
be desirable to have even more elementary approximations. For this purpose, quadratic approximations can be ex- 
ploited to further simplify equation (11). In particular, a candidate approximation is 

40 

m(t) = X{t- E[S 9 (t)))E[S(t)}^Var[S e (t)ElS(t)], (15) 

where 

45 

fI5 e (0]=qS(0 2 ]/2E[S(01 (16) 

and 

so 

Var(S c (t)) = (£[S(r) 3 J/3EIS(0]) - (E[S(t) 2 ]/2E[S(t)]) 2 % (17) 

with X*(t) being computed based on a suitably smoothed estimate of X. A convenient approximation for m(t) t which 
$ 5 tends to bo appropriate when X(t) changes slowly, is equation (15) without the second term, i.e., equation (15) with X" 
(t) = 0. See, S.G. Eick, W.A. Massey and W. Whitt, "The Physics of the M/G/~ Queue," Oper. fles., vol. 41 , pp. 731-742 
(1993). 

\ 
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An even simpler approximation is the infinite server system pointwise stationary approximation (PSA), namely 

m(t) = \(t)E{S(t)l (18) 

5 

The quadratic approximation in equation (15) differs from PSA in (18) by a time lag EJS^Ol and a space shift \*(t) Var 
[ 5 e(0]£lS(/)]/2. Equation (15) without the space shift is called shifted PSA. PSA in equation (18) is advantageously 
regarded as the natural initial approximation to m(t). The goal is to go beyond equation (18) and indicate when going 
beyond equation (18) is desirable. The refinement of equation (18) provided by equation (15) is only significantly better 

10 than (18) if the time lag and space shift are nonnegligible. Since the time lag and space shift can easily be computed 
explicitly, equation (15) is readily applied to determine when refinements to PSA should make a significant difference. 

For example, consider an operator staffing problem in which staffing changes are made every fifteen minutes. If 
the service times are approximately exponentially distributed with a mean of about one minute, then PSA should be 
adequate because E[S e (t] = E[S(t)] = 1 . However, if the average service time is fifteen minutes or more, as might occur 

15 with telephone calls handled by travel agents or emergency services (police, fire, etc.), then the refinements are likely 
to be more important when there is significant variation in X(t) over time. 

fn summary, equation (8) with the exact mean m(t) in equation(i) (and the variance as discussed below) is the 
infinite server approximation for the original finite-server problem. If a fixed exponential service-time distribution is 
assumed, then m(t) can be obtained for any arrival-rate function by solving the ODE in equation (14). If a phase-type 

20 service-time distribution is used, then m(t)car\ be obtained by solving a corresponding system of ODEs. Equation (15) 
is proposed as a convenient approximation to equation (1), and PSA in equation (18) as a convenient approximation 
to equation (15). The time lag and space shift in equation (15) help evaluate the importance of doing more than PSA. 

8 Considerations in determining v(t) 

25 

Recall that in step 330 in FIG. 3, an expression for the variance of the random variable QJt) was given in equations 
(2)-(4). The derivation of equations (2-4) are explained below, which derivation uses results from the analysis of sta- 
tionary infinite server systems. Recall that for a stationary infinite server system the peakedness is the ratio of the 
variance to the mean of the steady-state number of busy servers. Hence, without time dependence, and with the true 
30 peakedness, equation (2) would be exact, not approximate, by definition. A detailed, but complicated, expression for 
the peakedness is available in considerable generality, e.g., A.E. Eckberg, "Generalized Peakedness of Teletraffic 
Processes, " Proc. Tenth Int. Teletraffic Congress, Montreal, p. 4.4b. 3, June, 1983. Equation (3) is given as a useful 
approximation. 

There is additional support for the normal distribution with mean m(t) in equation (1 ) and variance v(t) in equation 
35 (2). For the special case of service-time distributions with a finite number of values, limiting results for stationary infinite 
server queueing systems extend to time-dependent arrivals to show that (t) is asymptotically normally distributed 
with (exact) mean m(t) in (6) and (exact) variance v(t) in equation (2) as the entire arrival rate function increases. See, 
P. W. Glynn and W. Whitt, "A New View of the Heavy-Traffic Limit Theorem for Infinite Server Queues, - Adv. Appl Prob 
vol. 23, pp. 188-209(1991). 

40 The simple formula for z(t) in equation (3) is obtained by assuming c\(t) and G/x)can be regarded as approximately 

stationary at time t. As for the approximation of the mean m(t) in equation (11), both cf (t) and GJx) are assumed to 
change much more slowly than X(t). Indeed, it is often advantageous to consider c?(t) constant. It is easy to obtain 
specific processes with c| (Yjconstant by considering a deterministic time change ofa stationary process with a well 
defined asymptotic variability parameter c a 2 (t). See, W.A. Massey and W. Whitt, "Unstable Asymptotics for Nonsta- 

4S tionary Queues," Math. Oper Res., vol. 19, pp. 267-291 (1994). 

As in the stationary case, z(t) = fcfft>+1)/2 when the service-time cumulative distribution function G/x) is expo- 
nential, and z(t)=c\(t) when the service-time cumulative distribution function is deterministic. In both cases, z(t)=*\ when 
c^ftfcl, as occurs with a Poisson process. 

The idea in (4) is to obtain c*(t) as the variance to mean ratio of the counting process in the neighborhood of t. 

so Look backward from f, since Q(t) depends on the arrivals before t. The estimate should be robust in the interval length 
x\ if A(t)-A(t<r\) is not too small and the variability is changing relatively slowly. If the arrival process is approximately a 
renewal process before time /, where interarrival times have squared coefficient of variation (SCV, variance divided by 
the square of the mean) c* then c\(t) should be just cf. 

The representation in equations (2)-(4) is convenient because modellers often have ideas about the relevant peak- 

55 edness z(t)or the asymptotic variability parameter c?(t) from extensive experience with stationary models. In that case, 
to directly calculate or estimate c\([) via equation (4) is straightforward. Thus, equations (2)-(4) permit the degree of 
stochastic variability to be specified directly in a relatively familiar way. 
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C. Considerations in refining s(t) 

Recall that in step 360 the inventive method optionally refines values for s(t). These refinements account for the 
fact that there are actually only finitely many servers. The goal is for equations (5a) and (5b) to hold for the random 
5 variable O c(t) (t) rather than the variable OJt), where Q s(t) (t) is the actual number of customers in the system (in service 
or waiting in the queue) when the number of servers as a function of time is given by the function s(t). 

The first refinement exploits a relationship between the tail probabilities P(Q„ >s) and P(Q^>s) in stationary systems 
with oo and s servers, respectively, when the offered load >£[S]\s large, where X is the arrival rate and E/S/is the mean 
service time. See, S. Halfin and W. Whitt, "Heavy-Traffic Limits for Queues with Many Exponential Servers." Oper. 
w Res. vol. 29, pp. 567-587 (1981); W. Whitt, "Understanding the Efficiency of Multi-Server Service Systems," Manage- 
ment Sci., vol. 38, pp. 708-723, (1992). In this regime, when PfQjzsj^a, then 

P(Q $ (t))>s = p 0 W«[1 + J2n z a (1 -a)exp(z^ /2)f 1 (19) 

15 

Hence when a target probability a is used in the basic infinite server system described above, the actual delay probability 
tends to be approximately p D (a.), which is somewhat greater than a, especially when a is not small. For example, when 
o=0.1 and 0.01, p D (a)=0. 132 and 0.115, respectively. 

Since p D (a) is a strictly increasing function of a, it is easy to construct the inverse function p D ~ 1 ®)\ i.e. if a= p D ' 1 
20 ft), then p D (a)=$. The refinement thus uses the infinite server system as described above with target probability p D ' 1 
(a) instead of a. 

The second refinement is an adjustment to s(t) after the function s(t) has been determined by the infinite server 
procedure, optionally with the first refinement above. The idea is to calculate PfQ^sft)) given that the function s(t) 
has been specified. The function s(t) is then adjusted, as necessary, so that, 

25 

PfQ^sft)) < a (20a) 



P(Q s(tH >s(t)-1)>a (20b) 

paralleling equations (5a) and (5b). For this purpose, it is advantageous to require that s(t) change only finitoly often. 
Then the computation is done at a representative set of time points t 

The procedure is iterative. Moreover, it is advantageous to proceed from earlier times to later times in making 
these adjustments because adjustments of $(t) at any time r can change the delay probabilities at later times. The 
procedure is repeated after a new function s(t) is determined. The iterative procedure stops when the function s(t) 
ceases to change significantly. 

If the arrival process is treated as Poisson and if the service times are treated as exponential or phase type, then 
the exact probabilities P(Q c(t £s(t)) required for the procedure above can be computed by solving a system of ODEs. 
(Exact computation is now possible because s(t) has been specified.) 

It is also advantageous to use another approximation called the modified-offe red-load (MOL) approximation for 
this purpose. The MOL approximation approximates P(Q s(t £k) for any k by the steady-state probability P(Q>k) as- 
suming that the number of servers is fixed at s, but uses m(t)/s(t) as the traffic intensity p, where m(t) is the mean 
number of busy servers in the infinite server system previously determined. The MOL approximation is easy to apply 
because the steady-state tail probability has a known explicit formula in terms of p. The MOL approximation is motivated 
by the fact that it is exact when the model is actually stationary. 

The MOL approximation for loss models is discussed in W. A. Massey and W. Whitt "An Analysis of the Modified 
Offered Load Approximation for the Nonstationary Erlang Loss Model," Ann. Appi Prob. vol. 4, pp. 1145-1160 (1994). 
The MOL approximation for delay models has essentially the same form (same idea works there). 

///. Conclusion 

The above method of determining the number of servers in a queueing system is not limited to the specific examples 
illustrated above. For example, the target probability a may be a function of time, in which case the extension is straight- 
forward in equation (8) where a is replaced by aft). Further, although the inventive method allows server staffing levels 
to be changed at arbitrary times, it is also possible to permit changes to made only at user specified times, t h /=1 : 2,... 
Then, throughout the interval ft, f /+ ,/ s(t) is set equal to the maximum value of s(t) as computed above over that interval. 
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The method also permits other constraints, such as servers being brought in or removed only in groups of certain sizes 
as by selecting s(t)\o be the next permissible group size In such cases. s(t) is the least integer such that all constraints 
are satisfied. Note that the inventive method is applicable when there is no waiting room or buffer space in the queue 
or, more generally, a finite waiting room, as well as when there is an infinite queue. Finally, the inventive method also 
applies in situations in which a serving system comprises either multiple queues (as in a supermarket) as well as when 
there is a single queue (as in many banks). 



Claims 



W 



15 



A method for operating a nonstationary finite server queuemg system, said method comprising the step of: 

operating said nonstationary finite server queueing system, wherein the number of servers in said nonsta- 
tionary finites server queueing system is determined by: 



a) determining a probability distribution function of a number of busy servers in a nonstationary infinite server 
queueing system as a function of time based on input parameters, wherein one of said input parameters 
characterizes an arrival rate of customers to said finite server queueing system, wherein one of said input 
parameters characterizes a service time distribution for a customer arriving at said finite server queueing 
system, wherein said probability distribution is characterized by a mean function and a variance function, 

20 wherein said mean function and said variance function are based on said input parameters, and 

b) determining a number of servers for said finite queueing system based on said distribution and based on 
a target probability. 



2. 



The method of claim 1 wherein said target probability is a delay probability characterizing the probability that a 
customer in said finite server queueing system will have to wait before receiving service. 

3. The method of claim 1 wherein said mean function is function of time, m(t), determined according to 

t 

«<0- fG c m (t-u)k(u)du. 

where G°Jt) m 1 -Gu(t) = P(S(u) > t) with S(u) being the service time of a customer arriving at time u to said finite 
server queueing system and X(u) is the arrival rate of said finite server queueing system at time u. 

4. The method of claim 3 wherein said variance function is a function of time, v(t), determined according to 

40 v(t) = z(t)m(t), 

where 



25 



30 



35 



45 



SO 



55 



£[5(0] 



gl Var{A(t)-A(t-T\)] 

[ k(u)du 
4-n 



where x\ > 0, where Aft) counts the number of arrivals before time t and where W?ris the variance. 



BNSOOCIO: <EP 0782066A2J. 



10 



EP 0 782 086 A2 

5. The method of claim 1 or claim 2 or claim 3 wherein said service time distribution is an exponential function char- 
acterized by a time-varying mean. 

6. The method of claim 1 or claim 3 wherein said arrival rate specifies a Poisson arrival process and wherein said 
5 variance function is equal to said mean function. 

7. The method of claim 1 wherein said probability distribution function of said number of busy servers in said non- 
stationary infinite server queueing system characterized as a normal distribution. 

w 8. The method of claim 1 wherein said number of servers for said finite server queueing system is adapted using a 
modified offered load approximation. 

9. The method of claim 1 wherein said target probability, a, is adapted so that the target probability is replaced by a 
quantity, p D 'Ha), where p D ' YaJ is the inverse of p D (a) where 

15 

p D (a) = [l . JTn z a (1-a)cxp(z^/2)] 1 

and where z w for a normal distribution characterized by a mean of zero and a standard deviation of one, N(0, 1), 
20 satisfies P(N(0, 1)> zj -a, where P(*) is the probability. 
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